Impacts of climate change on species distribution patterns of Polyspora sweet in China

Abstract Climate change is an important driver of species distribution and biodiversity. Understanding the response of plants to climate change is helpful to understand species differentiation and formulate conservation strategies. The genus Polyspora (Theaceae) has an ancient origin and is widely distributed in subtropical evergreen broad‐leaved forests. Studies on the impacts of climate change on species geographical distribution of Chinese Polyspora can provide an important reference for exploring the responses of plant groups in subtropical evergreen broad‐leaved forests with geological events and climate change in China. Based on the environmental variables, distribution records, and chloroplast genomes, we modeled the potential distribution of Chinese Polyspora in the Last Glacial Maximum, middle Holocene, current, and future by using MaxEnt‐ArcGIS model and molecular phylogenetic method. The changes in the species distribution area, centroid shift, and ecological niche in each periods were analyzed to speculate the response modes of Chinese Polyspora to climate change in different periods. The most important environmental factor affecting the distribution of Polyspora was the precipitation of the driest month, ranging from 13 to 25 mm for the highly suitable habitats. At present, highly suitable distribution areas of Polyspora were mainly located in the south of 25°N, and had species‐specificity. The main glacial refugia of the Chinese Polyspora might be located in the Ailao, Gaoligong, and Dawei Mountains of Yunnan Province. Jinping County, Pingbian County, and the Maguan County at the border of China and Vietnam might be the species differentiation center of the Chinese Polyspora. Moderate climate warming in the future would be beneficial to the survival of P. axillaris, P. chrysandra, and P. speciosa. However, climate warming under different shared socio‐economic pathways would reduce the suitable habitats of P. hainanensis and P. longicarpa.

With global warming, climate change will surpass habitat destruction and become the biggest threat to biodiversity in the future, which will have varying degrees of impact on biological individuals, populations, communities, and ecosystems (Bellard et al., 2012). According to the sixth assessment report of the Intergovernmental Panel on Climate Change (IPCC), the global average surface temperature in 2011-2020 was 1.09°C warmer than in 1850-1900; Over the next 20 years, global warming is expected to reach or exceed 1.5°C; By the end of the century, the global mean surface temperature will increase by 1.5-5.1°C (IPCC, 2021).
Existing studies have shown that climate warming will reduce the suitable area of some plants and show a fragmented distribution (Meng et al., 2019;Qiu et al., 2020;Wan et al., 2021). In addition, climate change may also have a positive impact on biodiversity.
Milder temperature and increased carbon dioxide may be beneficial to many plants, leading to accelerated biomass production.
Warmer winters may increase the survival rate of many species, and increased precipitation will benefit some moisture-loving plants (Bellard et al., 2012). Therefore, understanding the impact of future climate change on the distribution patterns of different species can provide a scientific basis for the formulation of species conservation strategies.
Niche model, also known as species distribution model, is a method of using known species distribution data and relevant environmental variables to build a model according to certain algorithm rules, judge the ecological needs of species, and project the calculation results to different time and space to predict the potential distribution area of species (Zhu et al., 2013). Neural Networks (ANN), and other models (Xu et al., 2015). Among them, MaxEnt model has short operation time, stable operation results  and high consistency between the predicted coverage and the actual distribution range of species (Booth et al., 2014;Ma et al., 2019). MaxEnt model is the most widely used niche model with good prediction effect (Giovanelli et al., 2008;Phillips & Dudik, 2008;Wang et al., 2020), which has been widely used in species evolution history (Jiang et al., 2016;Wang et al., 2021), impact of climate change on species distribution , endangered species protection (Wu et al., 2021), plant introduction and cultivation , and other fields.
Polyspora sweet belongs to Theaceae, species within the genus are evergreen tree or shrub. There are nearly 50 Polyspora species in the world, distributed in South and Southeast Asia, mainly in Malaysia, Indonesia, China, Vietnam, and other countries (Choo et al., 2020;Nguyet et al., 2020). Polyspora mainly grows in tropical and subtropical evergreen broad-leaved forests. There are six species of Polyspora in China, including P. axillaris, P. chrysandra, P. speciosa, P. hainanensis, P. longicarpa, and P. tiantangensis (Figure 1), which are mainly distributed in southwest and south China (Ming & Bartholomew, 2007), growing in mountain forests or shrubland, ridges, valleys, and mountainsides. The genus Polyspora has highornamental value, beautiful tree shape and blooming in winter.
It can be used as shade tree and street tree in gardens (Fan, Han, et al., 2021;Fan, Qian, et al., 2021;Ma et al., 2015). Some species have edible and medicinal value, the fruits contain natural antioxidants (Li et al., , 2019, and the extracts of roots and stems have cytotoxic activity (Fu, 2012;Tang, 2013;Xu et al., 2019).
Previous studies have shown that environmental and evolutionary factors play important roles in shaping species richness patterns of Theaceae in China (Rao et al., 2018). Polyspora has an ancient origin and completed species differentiation in the late Pliocene , and experienced the whole Quaternary climate change process. It is widely distributed and covers most of the subtropical areas in China, and is a typical representative in the subtropical evergreen broad-leaved forest. Therefore, the study on the distribution dynamics of Polyspora in different historical periods can not only provide new clues for understanding the species evolution of subtropical evergreen broad-leaved forest in China during the Quaternary ice age, but also provide scientific basis for the formulation of conservation measures and introduction and cultivation of Polyspora at present and in the future.

| Species distribution data
The geographical distribution data of Polyspora in China were mainly from the National Specimen Information Infrastructure (NSII, http://www.nsii.org.cn), Chinese Virtual Herbarium (CVH, https://www.cvh.ac.cn), Global Biodiversity Information Facility (GBIF, https://www.gbif.org), the published relevant literature, and field investigation from 2020 to 2022. For the sample points with a long history and no longitude or latitude information, Google Earth was used to locate and supplement them. Each specimen and each occurrence locality were carefully checked.
The misidentified specimen and occurrences recorded outside the native range of the species in the GBIF were deleted. To reduce the spatial self-correlation of population distribution points, avoid excessive fitting during MaxEnt operation, only one distribution point data was kept in the grid of 2.5 "× 2.5" on the map. Finally, in total 184 Chinese Polyspora occurrence sites were employed to build the model, including 63 P. axillaris, 43 P. chrysandra, 10 P. hainanensis, 19 P. longicarpa, 48 P. speciosa, and one P. tiantangensis (Table S1, Figure 2), these distribution records had covered the distribution range of each species. All the above, 184 distribution points were used in the prediction of the genus distribution area. According to unpublished data of our team, P. tiantangensis is actually an intraspecific variation type of P. longicarpa. In the simulation of species differentiation within the genus and the prediction at the species level, P. tiantangensis with only one distribution point was incorporated into P. longicarpa.

| Environmental data
Environmental factors such as climate, ultraviolet radiation, soil and terrain affect the growth, development and reproduction of plants, and then determine their distribution. In this study, five types of environmental datasets were selected, with a total of 36 environmental variables, including 19 bioclimatic variables, six UV-B variables, seven soil quality variables, three topographic variables, and one vegetation variable (Table S2).
Nineteen bioclimatic variables were downloaded from the World Climate Database (https://www.world clim.org). The WorldClim 1.4 dataset (Hijmans et al., 2005) was selected for paleoclimatic data, including the LGM and the mid-Holocene (MH; ~6000 years BP), the data were based on the Coupled Model Intercomparison Project Phase 5 (CMIP5). Since only CCSM4, MIROC-ESM and MPI-ESM-P Global Climate Models (GCMs) were available for the data of the LGM, to compare the prediction results of different GCMs, we selected the above three GCMs to simulate the ancient distribution of Polyspora. CCSM4 (The Community Climate System Model version 4) (Gent et al., 2011) is one of the most effective GCMs for predicting the impact of climate change on the distribution of animal and plant (Geng et al., 2022), and has the best precipitation prediction performance (Yang et al., 2020), especially for the precipitation prediction in southwest China (Yang, Yong, et al., 2021) (Watanabe et al., 2011), and a good prediction of rainfall in the Yangtze River basin of China (Yang et al., 2020). MPI-

ESM-P (Max Planck Institute for Meteorology Paleoclimate Model)
is a model specially designed for paleoclimate simulation (Braconnot et al., 2012;Jungclaus et al., 2013), which is accurate in simulating the trend of extreme temperature change in China (Jiang et al., 2017).
WorldClim 2.1 dataset (Fick & Hijmans, 2017) (Fan et al., 2020;, which improves the simulation capability of regional temperature and precipitation in China (Zhu et al., 2021). Four GCMs and four shared socio-economic pathways (SSPs) were used for future climate data. The four GCMs were the
CMCC-ESM2 and MPI-ESM1-2-HR are more accurate in predicting temperature in China, while EC-Earth3-Veg and CNRM-CM6-1 are more accurate in predicting precipitation in China (Yang, Zhou, et al., 2021;Zhu et al., 2021). There are four shared socio-economic pathways (SSPs) for future climate, including SSP126, SSP245, SSP370, and SSP585. SSPs can better reflect the correlation between socio-economic development and climate scenarios . Under these pathways, global warming will be 3-5°C by 2100 (Hausfather, 2020). The SSP126 scenario is a sustainable green path with global warming of 3-3.5°C by 2100; the SSP245 scenario is middle road, the world follows a path in which social, economic, and technological trends do not shift markedly from historical F I G U R E 2 Occurrence points of Polyspora in China. The black font in the map refers to the names of provinces, and the yellow fonts refers to the names of rivers.
patterns, with a warming range of 3.8-4.2°C by the end of the century; the SSP370 scenario is a regional competition route, with countries focusing on regional development, and a warming of 3.9-4.6°C by 2100; the SSP585 scenario envisages the world in which both economic output and energy consumption grow rapidly and without limit, with global warming of 4.7-5.1°C by 2100.
In total, six UV-B radiation variables were obtained from the Global UV-B radiation database (gIUV, https://www.ufz.de/gluv/) (Beckmann et al., 2014). Seven soil quality data were downloaded from the Harmonized World Soil Database v 1.2 (HWSD, https:// www.fao.org/soils -porta l/data-hub/soil-maps-and-datab ases/ harmo nized -world -soil-datab ase-v12/zh/) (Fischer et al., 2008).  (Shabani et al., 2019;Wu et al., 2015Wu et al., , 2021. The temperature value in the original layer of paleoclimatic variables were reduced by 10 times with the R package "terra" (Robert, 2022), which was unified with the current and future climate layers. All the environmental variables were trimmed to China region using R package "terra" (Robert, 2022), resamped to a unified spatial resolution of 2.5′, and converted to ASC format. China region was based on the Chinese vegetation layer obtained by RESDC. All the coordinate systems used in this study were the World Geodetic Coordinate System 1984 (WGS1984).
To avoid over-fitting or inaccurate modeling of MaxEnt model due to strong correlation between environmental variables (Hu & Liu, 2014;Li et al., 2020), the correlation between climate factors were comprehensively tested by using R package "corrplot" (Wei & Simko, 2021) Figure 4, Appendix A) based on the standard protocol of species distribution model (Zurell et al., 2020). The entire workflow of data source and preprocessing, modeling parameter generation, model calculation, and performance verification were described in detail.

| MaxEnt species distribution modeling
Geographical distribution data of Polyspora species and the climate data of four periods were imported into MaxEnt v3.4.4, 75% of the distribution points were randomly selected for modeling, and the remaining 25% were used for verification. The maximum number of iterations was set to 1000, background points 10,000, replications bootstrap 10, and use the average value of 10 operations to ensure the reliability of the results. The R package "ENMeval 2.0" (Kass et al., 2021) was used to optimize two key modeling parameters: regularization multiplier (RM) and feature combination (FC). For Polyspora and P. axillaris, RM was set to 0.5-4, with an interval of 0.5 each time, a total of 8 values. The RM of P. chrysandra, P. Hainanensis, P. longicarpa, and P. speciosa were set to 1-4 with each interval of 0.5, a total of 7 values. RM = 0.5 had been tested for these species, but the results were not desirable, and the difference between the running results of 0.5 and 1 was too large. We selected six types for FC testing, which were L, LQ, LQH, H, LQHP, and LQHPT. Among them, L means linear, Q means quadratic, H means hinge, P means product, and T means threshold. Therefore, a total of 42 to 48 sets of parameter combinations were used for the ENMeval test. For Polyspora, P. axillaris, P. chrysandra, and P. speciosa, the "block" partition method was used to perform ENMeval operation to limit the autocorrelation in the large scale space. Since the sample sizes of P. hainanensis and P. longicarpa were <25, the "Jackknife" method was used for ENMeval operation to generate the maximum available information model and improve the prediction accuracy (Galante et al., 2018;Pearson et al., 2007;Shcheglovitova & Anderson, 2013).
The model performance was evaluated according to akaike information criterion (AICc), delta akaike information criterion (delta AICc), the difference between training areas under curves (AUCs) and testing AUC (AUC.diff), and the 10% training commission rate (OR10).
The model with the lowest AICc (delta.AICC = 0) and the higher AUC value (AUC > 0.9) was considered to be the best model. Since models with delta.AICc < 2 are all reliable (Phillips et al., 2017), to improve the reliability of models, we also considered models with smaller AUC. Diff and OR10 for some species (Warren & Seifert, 2011).
The simulation accuracy was evaluated by the value of the area under the receiver operating characteristic curve (AUC) (Lobo et al., 2008) and the continuous Boyce index (CBI) (Boyce et al., 2002).
The range of AUC values is [0, 1], and the closer it is to 1, the more accurate the simulation results are (Jiang et al., 2016). AUC > 0.9 indicates that the simulation result is very accurate, and the simulation results with AUC < 0.7 are not credible (Elith et al., 2006).

| Phylogenetic analysis and differentiation time estimation
To better reveal the phylogeny and species differentiation relationship of Chinese Polyspora, and confirm the simulation results of MaxEnt model, the complete chloroplast genomes of six Polyspora species were downloaded from GenBank, of which two species (P. tiantangensis and P. chrysandra) were assembled, annotated and uploaded by our team, the chloroplast genome sequence of Apterosperma oblata was also downloaded for used as an outgroup in the phylogenetic tree (Table S4). We aligned the seven complete chloroplast genomes with MAFFT v7.450 (Katoh & Standley, 2013), selected the best nucleotide substitution model (GTR + I + G) with modelgenerator v0.85 (Keane et al., 2006), then constructed the   Table S5). The optimized parameters significantly reduced the model complexity and were more suitable for modeling migration in different periods. Therefore, RM = 1.5 and FC = LQH were selected as modeling parameters for niche simulation of Polyspora in this study. Similar methods were used to select the optimal modeling parameters for niche simulation of each species in the genus Polyspora. The results were shown in Tables S3 and   S6-S10, and Figures S6-S10.

| Optimal model and model accuracy evaluation
The selected parameter combinations were used to simulate and predict the historical, contemporary, and future distribution areas of Polyspora species. The AUC values of the training set and test set in each period, each GCM and each SSP were all >0.9, and the CBI values were all >0.5, indicating that the model had good-fitting effect and high-prediction accuracy.

| Potential suitable habitat of Polyspora under paleoclimate scenarios
The distribution simulation of Polyspora during LGM period using three GCMs showed that the highly suitable habitats were mainly located in Yunnan Province. In CCSM4 and MIROC-ESM, the highly suitable habitat was mainly located in southwestern Yunnan. While in MPI-ESM-P, the highly suitable area was significantly shifted to the eastern Yunnan (Figure 9).
In the MH, the habitats of Polyspora under the three GCMs were close to the current suitable habitats, with almost overlapping centroids ( Figure 10), but the suitable habitat area was less than current.
Especially in MPI-ESM-P, the highly suitable habitat area was only 29,854 km 2 (Figure 9, Table S12). Compared with the LGM, the highly suitable distribution area expanded eastward, and the distribution centroid also shifted eastward ( Figure 10).  Figure 11, Table S12).

| Potential suitable habitat of
Under the SSP585 scenario, except for the EC-Earth3-Veg, the highly suitable area predicted by the other three GCMs would be lost to varying degrees, with loss rates ranging from 12.95% (MPI-ESM1-2-HR) to 36.71% (CNRM-CM6-1) ( Figure 11,  (Figures 11 and 12). It could be seen from the centroid distribution map that under different SSPs in the future, the distribution centroids of Polyspora would all move southeast to Guangxi. The location of the centroids under different SSPs were very close ( Figure 10).

| Simulation of distribution patterns of Polyspora species in different periods
The current potential highly suitable areas of five species of Polyspora (P. axillaris, P. chrysandra, P. speciosa, P. hainanensis, and P. longicarpa) were species-specific ( Figure 13, Figure S11). Highly suitable areas of P. axillaris were located in southeastern Guangdong and Taiwan. The highly suitable area of P. hainanensis was only located in the Hainan Island. Highly suitable areas of P. chrysandra and P. longicarpa were both located in the west and southwest Yunnan, and the two species have a large area of overlapping distribution. Polyspora speciosa had the widest range of high suitability, which were located in southeast Sichuan, Chongqing, Guizhou, northern Guangxi and Taiwan.
Polyspora axillaris and P. speciosa overlapped slightly in Taiwan.
Since the LGM, P. axillaris had experienced a large expansion, and the optimal distribution area in the MH was slightly higher than that in the current period ( Figure S12). At present, highly suitable habitat for P. axillaris was located in the Pearl River Delta and Taiwan ( Figure S13). In the future, under different SSPs, the total suitable area and highly suitable area would both increase significantly, and expand eastward and westward at the same time. Under EC-Earth3-Veg_SSP585 scenario, the highly suitable area would reach 110,746 km 2 in 2100 (Figures S14 and S15, Table S13).
Under the three GCMs of paleoclimate, the potential suitable areas of P. chrysandra during the LGM were more than that of current period ( Figure S16, Table S14). From LGM to MH, the suitable habitat decreased, and then expanded after the MH. Currently, highly suitable areas were located in southwestern and southern Yunnan ( Figure S17). In the future, except for the CMCC-ESM2_SSP245 scenario, the potential suitable habitat of P. chrysandra would continue to expand, with an average increase of 40,000 km 2 by 2100, but the average loss of highly suitable habitat was 10,000 km 2 (Table S14, Figures S18 and S19).

F I G U R E 8 Potential distribution areas of Chinese
Polyspora under current environmental conditions. When presence probability is <0.1, unsuitable region; when presence probability is 0.1-0.3, lowly suitable region; when presence probability is 0.3-0.5, moderately suitable region; and when presence probability is >0.5, highly suitable region.
Highly suitable areas of P. speciosa during the LGM were mainly located in southeast Sichuan, Chongqing and northern Guizhou, and expanded to Guangxi in the MH period ( Figure S20). At present, in addition to the above areas, there were also highly suitable habitats for P. speciosa in central and eastern Taiwan, and a large range of suitable habitats in eastern and southern China ( Figure S21). In the future, under different SSPs, the total suitable area would continue to expand, but the highly suitable area would be shrunken in a small range under the CMCC-ESM2 and CNRM-CM6-1 modes ( Figures S22-S25, Table S15).
Current suitable habitats of P. hainanensis were mainly located in southern Hainan and eastern Taiwan, and the total suitable areas of P. hainanensis were less than other Chinese Polyspora species, which were 37,158 km 2 ( Figure S26, Table S16). The simulation results of three GCMs of paleoclimate showed that during the LGM, there was no moderately and highly suitable habitat for P. hainanensis, only lowly suitable habitat was found in southern Yunnan ( Figure S27).
In the MH, highly suitable habitat of P. hainanensis appeared in southern Hainan, while the suitable habitat in Yunnan was lost. With the future climate warming, the total suitable area of P. hainanensis would be reduced, but the highly suitable area would be relatively stable ( Table S16, Figures S28 and S29).
During the LGM, highly suitable habitats of P. longicarpa were located in most areas of Yunnan and southern Taiwan, and under the LGM_MPI-ESM-P model, there were also highly suitable habitats in southeastern Tibet ( Figure S30). In the MH, the total suitable area decreased, but the highly suitable areas in the two paleoclimatic periods were more than that in the present period ( Figure S30, Table S17). At present, P. longicarpa was mainly distributed in the western Yunnan and central Taiwan, with highly suitable area of 72,858 km 2 (Table S17, Figure S31). Under various GCMs and SSPs combination in the future, except for the MPI-ESM1-2-HR_SSP126 model, the total and highly suitable areas of P. longicarpa would shrink westward by 2100. Under the CNRM-CM6-1_SSP585 model, only 18,536 km 2 would remain in the highly suitable area, with a loss rate of 74.56% ( Figures S32 and S33, Table S17).

F I G U R E 9
Potential distribution areas of Chinese Polyspora under paleoclimate scenarios. Comparison of three global climate models (CCSM4, MIROC-ESM, and MPI-ESM-P) in two paleoclimatic periods, the last glacial maximum and mid-Holocene.

| Phylogenetic analysis of the Chinese Polyspora based on the chloroplast genomes
The phylogenetic tree based on chloroplast genomes showed that the earliest branch of Chinese Polyspora was P. speciosa, and then it diverged into two big branches. Divergence time by Beast showed that the divergence between P. axillaris and P. hainanensis occurred at 3.12 Ma. At 5.58 Ma, P. chrysandra differentiated with P. longicarpa and P. tiantangensis. The differentiation of P. longicarpa and P. tiantangensis was the latest, and the divergence occurred during the Last Interglacial (LIG), 164 ka bp ( Figure 14).

| Key environmental factors shaping species distribution
Factors affecting species distribution include environmental factors, biological factors, and the characteristics of species themselves . Environmental factors mainly affect species distribution at large spatial scales (King et al., 2020). Biological factors and the characteristics of species themselves mainly affect the distribution of species in a relatively small scale (Zhu et al., 2013). The results of this study showed that Polyspora were widely distributed in the south of the Yangtze River, and the highly suitable area was mainly located in the south of 25° north latitude (Figure 8, current). The climate type of this area is mainly subtropical monsoon humid climate, and the vegetation type is mainly subtropical evergreen broad-leaved forest, Theaceae is one of the typical representative plant group. in the distribution of P. axillaris was annual precipitation, the annual precipitation in the Pearl River Delta region was 1600-2300 mm, and the annual precipitation in Taiwan Island was over 2500 mm, which created the best conditions for the growth of P. axillaris. UV-B radiation has a significant impact on aboveground organs of plants, thereby limiting the distribution of species (Wu et al., 2021). In addition to bioclimatic factors, UV-B radiation also played an important role in the distribution of P. axillaris. Isothermality was the primary ecological factor determining the distribution of P. chrysandra and P. longicarpa, these two species were mainly distributed in southern and western Yunnan, where the annual temperature difference was small and the monthly temperature varied greatly, providing a high isothermality for the species distribution. In addition to bioclimatic factors, eleva-

| Species differentiation and glacial refugia speculation
In the evolutionary history of species, niche differentiation can promote species differentiation (Zhu et al., 2013), and species distribution models can test the hypotheses in evolutionary biology (Peterson & Nyari, 2008;Ree et al., 2005). Based on the phylogeny and divergence time of the chloroplast genome, species differentiation of genus Polyspora began in the late Miocene, and most species completed their differentiation in the late Pliocene ( Figure 14).
Our estimated differentiation time was consistent with the research results of Yu et al. (2017). Since the Pliocene, the surface temperature had generally declined ( Figure 14). During the LGM, some differentiated species stayed in situ or migrated to the common refugia. to Vietnam, were the common suitable distribution areas for four species of Polyspora (P. chrysandra, P. speciosa, P. hainanensis, and P. longicarpa), which were located in four clades of the phylogenetic tree ( Figure 14). We speculated that this area might be the differentiation center of Polyspora species in China. The LGM was the closest to the human environment in the last 20,000 years, and had a great contrast with the modern climate. At that time, the climate in China was cold and dry, and most areas were grassland and desert. The edge of the grassland could reach the north of the modern evergreen broad-leaved forest, while the evergreen broad-leaved forest reached the modern tropical area, and the tropical forest completely disappeared . Figures S34-S38 shows the niche overlaps between the LGM and the current F I G U R E 1 3 Potential highly suitable habitat of Polyspora species in China. Different colors represent highly suitable habitats for different species, among which yellow and orange indicate highly suitable habitats shared by two species.
F I G U R E 1 4 Phylogenetic relationships and divergent time estimated of Chinese Polyspora inferred from maximum likelihood analysis based on chloroplast genomes. The upper part shows the mean temperature difference between ancient and present. The middle part is the phylogenetic tree based on the maximum likelihood method, the numbers on the nodes indicate bootstrap values. The lower part is the divergent time of Chinese Polyspora, dark blue bars represent the node ages of 95% confidence interval.
period of Polyspora species in China. During the LGM, suitable habitat of P. axillaris was mainly located in the Lianhua Mountain area in eastern Guangdong, where might be the glacial refugium of P. axillaris ( Figure S34). Polyspora chrysandra was still widely distributed during the LGM, and the glacial refugium might be located in Ailao Mountain ( Figure S35). The glacial refugium of P. speciosa might have been located around Chongqing, and expanded to southeast after the glaciation ( Figure S36). Due to the disappearance of tropical forests during the LGM, P. hainanensis lost its habitat and migrated northwest. Dawei Mountain in southeastern Yunnan might be its glacial refugium ( Figure S37). The glacial refugium of P. longicarpa might be located in the Gaoligong Mountain of western Yunnan, where it expanded northward in a small range after the glacial period ( Figure S38).
The unique topography and climate in southwest China provide rich habitats for animals and plants, making it not only a hotspot of biodiversity in the world, but also an ideal place of species origin and glacial refugia Qiu et al., 2011;Tang et al., 2018).
By comparing the distribution ranges of various species of Polyspora in current and LGM, we found that the main glacial refugia of the genus were located in southwestern China. In addition, Lianhua Mountain in Guangdong Province also formed a relatively independent niche, which weakened the impact of Quaternary climatic turbulence to a certain extent, and became another potential suitable area of P. axillaris during the LGM.

| Impacts of future climate change on the distribution of Polyspora
Climate change is the major factor affecting the large-scale distribution pattern and migration pattern of species (Pearson & Dawson, 2003;Wu et al., 2021). The prediction results of several climate models show that the surface temperature in China will continue to rise in the future, and the annual precipitation in most parts of the country will also increase. By the end of the century, the average temperature will increase by 2.7-5.4°C, and the annual average precipitation will increase by 17%-30% (Yang, Zhou, et al., 2021). The prediction results showed that the total suitable area of Polyspora in China would generally decrease under the four trend. In general, some Polyspora species will lose part of their suitable area, but each species has a relatively stable habitat area. Climate change will not cause large-scale migration or extinction of Polyspora species, but reasonable control of carbon emissions will beneficial to the survival and distribution of Polyspora species.

| CON CLUS ION
In this study, the species geographical distribution data and environmental variables were used to simulate the distribution dynamics of

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in Dryad at https://doi.org/10.5061/dryad.jq2bv q89x.  Barstow, M., & Rivers, M. (2017). The Red List of Theaceae.

Model objective
Model objective: Forecast and transfer.
Target output: Habitat suitability maps, continuous occurrence probabilities, and binary maps of potential presence.

Location
Location: China. Temporal resolution: Averages over each period.

Scale of analysis
Boundary: political.

Biodiversity data
Observation type: field survey, GPS tracking, citizen science.
Response data type: presence/absence.

Hypotheses
Hypotheses: Effects of environmental variables (bioclimatic, light, soil, topographic, vegetation) on species distribution and differentiation.

Assumptions
Model assumptions: Sampling is adequate, and representative, spatial sampling bias has been successfully corrected. The distribution and differentiation of Polyspora were driven by environmental variables. Polyspora species are at equilibrium with their environment, their niches expand, contract, or migrate with climate change.

Algorithms
Modelling techniques: MaxEnt. Clipping: We clipped all data to the political boundary of China.
Scaling: Duplicate records in the same localities were removed and spatial autocorrelation was minimized by randomly removing occurrences within 5 km of each other.
Cleaning: We carefully verified each specimen and removed the misidentified occurrence localities.
Absence data: Not applicable.
Background data: We generated 10,000 random background points within the study area.
Errors and biases: We checked each occurrence location coordinate one by one, removed error values, outliers, duplicates, and insufficient information items, so the error rate deemed low.

Data partitioning
Training data: We randomly selected 75% of the data for model building and 25% for validation of the predictions.
Validation data: Cross-validation method.
Test data: Not applicable. Spatial resolution: Data were collected and prepared at the same resolution as the bioclimatic data: 2.5 min spatial resolution (approximately 5 km).
Temporal resolution: The raw resolution of bioclimatic data was 2.5 min (approximately 5 km); the raw resolution of gIUV database was 5 arc-min (approximately 10 km); the raw resolution of HWSD database was 30 arc-s (approximately 1 km); the raw resolution of the vegetation data was 1 km.
Data processing: For UV-B, soil and vegetation variables, we used raster package (version 3.5-15) to unify the resolution consistent with bioclimatic variables.
Errors and biases: Not applicable.
Dimension reduction: Predictor variables were standardized and used correlation analysis to avoid highly correlated variables. We filtered 14-19 environment variables for modeling. Data processing: The paleoclimatic temperature grid layer (bio 1bio 2, bio 4-bio 11) was uniformly reduced by 10 times using the terra package (version 1.5-12) to be consistent with the data units (°C) of other periods.
Quantification of novelty: Distance to training data.

Variable pre-selection
Variable pre-selection: Ecological pre-selection of variables we deemed important for the species, down to 14-19 predictors.

Model estimates
Coefficients: We used the area under the curve (AUC) and continuous Boyce index (CBI) to evaluate model robustness and excluded models for which the AUC was below 0.9 and for which the CBI was below 0.5.
Parameter uncertainty: Cross-validation was used to determine the optimum tree size yielding the most robust predictions.
Variable importance: Jackknife method in the model.

Model selection-model averaging-ensembles
Model selection: Selection of the best individual models based on the AUC value (>0.9) and CBI value (>0.5).
Model averaging: Predicted mean of selected individual models.

Analysis and correction of non-independence
Spatial autocorrelation: None.
Nested data: None.

Threshold selection
Threshold selection: The suitability grades of Polyspora were divided into four regions: unsuitable region (0-0.1), lowly suitable region (0.1-0.3), moderately suitable region (0.3-0.5), and highly suitable region (0.5-1). In the simulation of centroid transfer, suitable area change, and interspecific overlapping distribution area, the commonly used fixed threshold of 0.5 were used as the critical value of species distribution/nondistribution.

Performance statistics
Performance on training data: AUC, BIC.
Performance on validation data: AUC, BIC.
Performance on test data: AUC, BIC.

Plausibility check
Response shapes: We checked model plausibility by assessing partial dependence plots.
Expert judgement: Maps of modelled predictions were checked by experts for an ad-hoc subset of species.

PR ED I C TI O N Prediction output
Prediction unit: Continuous habitat suitability index (0-1), gain and loss habitat estimation from the binary prediction (presence/ absence).

Uncertainty quantification
Algorithmic uncertainty: Not applicable.
Input data uncertainty: Not applicable.
Parameter uncertainty: Not applicable.
Scenario uncertainty: We compared the differences in species suitable habitats in different GCMs during the same period, focusing on the same distribution areas to reduce the uncertainty in scenarios.
Novel environments: Climate map visualization for each period and scenario.